home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / calcprf1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  3.1 KB  |  106 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include "clustalw.h"
  6.  
  7. /*
  8.  *   Prototypes
  9.  */
  10.  
  11. extern void *ckalloc(size_t bytes);
  12. extern void ckfree(void *);
  13.  
  14. void calc_prf1(int **profile, char **alignment, int *gaps,
  15.   int matrix[NUMRES][NUMRES],
  16.   int *seq_weight, int prf_length, int first_seq, int last_seq);
  17.  
  18. /*
  19.  *   Global variables
  20.  */
  21.  
  22. extern int max_aa,gap_pos1,gap_pos2;
  23.  
  24. void calc_prf1(int **profile, char **alignment, int *gaps,
  25.   int matrix[NUMRES][NUMRES],
  26.   int *seq_weight, int prf_length, int first_seq, int last_seq)
  27. {
  28.  
  29.   int **weighting; 
  30.   int f, sum2;        
  31.   float scale;
  32.  
  33.   int i, d, pos, res, count;
  34.   int   r, numseq;
  35.  
  36.   weighting = (int **) ckalloc( (NUMRES+2) * sizeof (int *) );
  37.   for (i=0;i<NUMRES+2;i++)
  38.     weighting[i] = (int *) ckalloc( (prf_length+2) * sizeof (int) );
  39.  
  40.   numseq = last_seq-first_seq;
  41.  
  42.   sum2 = 0.0;
  43.   for (i=first_seq; i<last_seq; i++)
  44.        sum2 += seq_weight[i];
  45.  
  46.   for (r=0; r<prf_length; r++)
  47.    {
  48.       for (d=0; d<=max_aa; d++)
  49.         {
  50.             weighting[d][r] = 0;
  51.             for (i=first_seq; i<last_seq; i++)
  52.                if (d == alignment[i][r]) weighting[d][r] += seq_weight[i];
  53.         }
  54.       weighting[gap_pos1][r] = 0;
  55.       for (i=first_seq; i<last_seq; i++)
  56.          if (gap_pos1 == alignment[i][r]) weighting[gap_pos1][r] += seq_weight[i];
  57.       weighting[gap_pos2][r] = 0;
  58.       for (i=first_seq; i<last_seq; i++)
  59.          if (gap_pos2 == alignment[i][r]) weighting[gap_pos2][r] += seq_weight[i];
  60.    }
  61.  
  62.   for (pos=0; pos< prf_length; pos++)
  63.     {
  64.       if (gaps[pos] == numseq)
  65.         {
  66.            for (res=0; res<=max_aa; res++)
  67.              {
  68.                 profile[pos+1][res] = matrix[res][gap_pos1];
  69.              }
  70.            profile[pos+1][gap_pos1] = matrix[gap_pos1][gap_pos1];
  71.            profile[pos+1][gap_pos2] = matrix[gap_pos2][gap_pos1];
  72.         }
  73.       else
  74.         {
  75.            scale = (float)(numseq-gaps[pos]) / (float)numseq;
  76.            for (res=0; res<=max_aa; res++)
  77.              {
  78.                 f = 0.0;
  79.                 for (d=0; d<=max_aa; d++)
  80.                      f += (weighting[d][pos] * matrix[d][res]);
  81.                 f += (weighting[gap_pos1][pos] * matrix[gap_pos1][res]);
  82.                 f += (weighting[gap_pos2][pos] * matrix[gap_pos2][res]);
  83.                 profile[pos+1][res] = (int  )(((float)f / (float)sum2)*scale);
  84.              }
  85.            f = 0.0;
  86.            for (d=0; d<=max_aa; d++)
  87.                 f += (weighting[d][pos] * matrix[d][gap_pos1]);
  88.            f += (weighting[gap_pos1][pos] * matrix[gap_pos1][gap_pos1]);
  89.            f += (weighting[gap_pos2][pos] * matrix[gap_pos2][gap_pos1]);
  90.            profile[pos+1][gap_pos1] = (int )(((float)f / (float)sum2)*scale);
  91.            f = 0.0;
  92.            for (d=0; d<=max_aa; d++)
  93.                 f += (weighting[d][pos] * matrix[d][gap_pos2]);
  94.            f += (weighting[gap_pos1][pos] * matrix[gap_pos1][gap_pos2]);
  95.            f += (weighting[gap_pos2][pos] * matrix[gap_pos2][gap_pos2]);
  96.            profile[pos+1][gap_pos2] = (int )(((float)f / (float)sum2)*scale);
  97.         }
  98.     }
  99.  
  100.   for (i=0;i<=max_aa;i++)
  101.     ckfree((void *)weighting[i]);
  102.   ckfree((void *)weighting);
  103.  
  104. }
  105.  
  106.